#ifndef _RANDOM_H_
#define _RANDOM_H_

namespace math 
{
	struct distribution_tag {};
	struct uniform_tag : public distribution_tag {} ;
	struct normal_tag : public distribution_tag {} ;
	struct laplace_tag : public distribution_tag {} ;

	template <typename value_type, typename dis_tag, int low_bound, int up_bound>
	struct random_generator ;
	

	template <int low_bound, int up_bound> struct random_generator <double, uniform_tag, low_bound, up_bound> 
	{
	public :
		typedef double value_type ;
		value_type seed (unsigned s) {srand (s);}
		value_type random() const { return  low_bound + value_type(rand()) * (up_bound - low_bound) / value_type (RAND_MAX); }
	};

	template <int low_bound, int up_bound> struct random_generator <float, uniform_tag, low_bound, up_bound> 
	{
	public :
		typedef float value_type ;
		value_type seed (unsigned s) {srand (s);}
		value_type random() const { return  low_bound + value_type(rand()) * (up_bound - low_bound) / value_type (RAND_MAX); }
	};

	template <int low_bound, int up_bound> struct random_generator <int, uniform_tag, low_bound, up_bound> 
	{
	public :
		typedef int value_type ;
		value_type seed (unsigned s)  {srand (s);}
		value_type random() const { return  low_bound + rand() % (up_bound - low_bound); }
	};

	template <int low_bound, int up_bound> struct random_generator <double, normal_tag, low_bound, up_bound> 
	{
		typedef double value_type ;
		value_type seed (unsigned s) {srand (s);}
		typedef random_generator <double, uniform_tag, 0, 1> uniform_generator_type ;
		uniform_generator_type uniform_generator ;
		random_generator (value_type _m, value_type _s): m(_m), s (_s) {}

		// mean m, standard deviation s 
        // normal random variate generator 
		value_type random()	const 
		{				        
			value_type x1, x2, w, y1;
			static value_type y2;
			static int use_last = 0;
			if (use_last)		        // use value from previous call 
			{
				y1 = y2;
				use_last = 0;
			}
			else
			{
				do {
					x1 = 2.0 * uniform_generator.random() - 1.0;
					x2 = 2.0 * uniform_generator.random() - 1.0;
					w = x1 * x1 + x2 * x2;
				} while ( w >= 1.0 );

				w = sqrt( (-2.0 * log( w ) ) / w );
				y1 = x1 * w;
				y2 = x2 * w;
				use_last = 1;
			}
			return( m + y1 * s );
		}
	private : 
		const value_type m , s; 
	};

}
#endif
